#################################
##' Code for:
##' Title: "Effective climate clubs require ambition, leverage, and insulation"
##' Author: Sam S. Rowan
##' Contact: sam.rowan [at] concordia.ca
##' Date: 2023-10-19
#################################
##' This script:
##' Illustrating accession
#################################

## Working directory
# setwd("/Users/srowa/Dropbox/Projects/climate_clubs_trade")
# -- Set working directory to replication folder


## Need objects from '05_ideal_type_clubs.R'
source("06_ideal_type_clubs.R")


#### Illustrating accession ####

# Create a dummy variable for "accession" to EU club

eu_accession_club <- total_clubs %>% 
  rename(reporter_in_club = reporter_in_eu, 
         partner_in_club = partner_in_eu) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>% 
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club) %>%
  distinct() %>% 
  as.data.frame() %>% 
  # Merge in GHG emissions data on reporter
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  filter(!is.na(reporter_in_club),
         !is.na(ghg_multiplicative_log)) %>% 
  # Create 'accede' variable
  mutate(reporter_accede = ifelse(ghg_multiplicative_log*.1 <= export_exposure_to_club, 
                                  1, 0)) %>% 
  mutate(reporter_in_accede = ifelse(reporter_in_club==1, 1, 
                                     ifelse(reporter_accede == 1, 1, 0))) %>% 
  select(iso3c = reporter_iso, accession_club = reporter_in_accede)

# Re-calculate trade exposure to EU and accession club
accession_eu_exposure <- total_clubs %>% 
  rename(reporter_in_club = reporter_in_eu, 
         partner_in_club = partner_in_eu) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>% 
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  select(reporter_iso, partner_iso, export_share_to_partner, exports_world,
         reporter_in_club, partner_in_club, export_exposure_to_club) %>% 
  # Merge in accession club on partner
  full_join(., eu_accession_club, by = c("partner_iso" = "iso3c")) %>% 
  rename(partner_in_accession = accession_club) %>% 
  full_join(., eu_accession_club, by = c("reporter_iso" = "iso3c")) %>% 
  rename(reporter_in_accession = accession_club) %>% 
  filter(!is.na(reporter_in_accession),
         !is.na(partner_in_accession)) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_accession) %>% 
  mutate(export_exposure_to_accession = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club,
         reporter_in_accession, export_exposure_to_accession) %>%
  distinct() %>% 
  as.data.frame() %>% 
  mutate(membership = ifelse(reporter_in_club == 1, "Core", 
                             ifelse(reporter_in_accession == 1 & reporter_in_club == 0, "Accede", 
                                    "Non-member")))

knitr::kable(table(accession_eu_exposure$membership))
# Recall, no data for Czechia EU-27 -> 'EU-26'

# Average shift in trade exposure
accession_eu_exposure %>% 
  group_by(membership) %>% 
  summarize(exposure_shift = mean(export_exposure_to_accession - export_exposure_to_club))


## Accession: scatterplot
p_accede1 <- accession_eu_exposure %>%
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  ggplot(., aes(x = ghg_multiplicative_log, y = export_exposure_to_accession,
                label = reporter_iso, color = membership)) +
  geom_smooth(method = "lm", se = F) +
  geom_smooth(aes(x = ghg_multiplicative_log, y = export_exposure_to_club),
              method = "lm", se = F, linetype = 2) +
  geom_text(size = 5) +
  geom_point(aes(x = ghg_multiplicative_log, y = export_exposure_to_club),
             size = 1, pch = 1) +
  geom_segment(aes(x = ghg_multiplicative_log, xend = ghg_multiplicative_log,
                   y = export_exposure_to_club, yend = export_exposure_to_accession),
               size = 0.75) +
  geom_abline(slope = 0.1, intercept = 0, linetype = 3) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  annotate("text", x = 10, y = 1.0, label = "Accession rule", adj = 1, size = 5) +
  # scale_y_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_color_manual(values = c("#79b321", "#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(y = "Export exposure to club members",
       color = "Club member",
       x = expression(paste("GHG emissions" %*% "GHG emissions per capita, log transformed")))

## Accession: boxplot
p_accede2 <- accession_eu_exposure %>%
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  ggplot(., aes(export_exposure_to_accession,
                group = membership,
                fill = membership)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(breaks = c(-0.25, 0, 0.25), labels = c("Accede", "Core", "No")) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  scale_fill_manual(values = c("#79b321", "#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(x = "",
       fill = "Club member",
       y = "Club member") 

cowplot::plot_grid(p_accede1, p_accede2, align = "v", ncol = 2, rel_widths = c(2/3, 1/3))
ggsave("text/comtrade_accede_stacked.pdf", units = "in", height = 6, width = 12)



## How many new countries would join after the first accession
accession_eu_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_accede_2 = ifelse(ghg_multiplicative_log*.1 <= export_exposure_to_accession, 
                                    1, 0)) %>% 
  mutate(true_accession_2 = ifelse(reporter_accede_2==1 & membership == "Non-member", 1, 0)) %>% 
  summarize(table(true_accession_2)) # One

